﻿using System;
using System.Collections.Generic;
using System.Text;
using System.IO;

namespace EquationSolver
{
    /// <summary>
    /// Represents a Gauss-Jordan linear equations solver
    /// </summary>
    public class EquationSolver
    {
        //Data members
        private List<double[]> _matrix = null; //Matrix of equations
        private int _numOfVars = -1; //Number of variables

        /// <summary>
        /// Constructor
        /// </summary>
        public EquationSolver()
        {
            _matrix = new List<double[]>();
        }

        /// <summary>
        /// Adds an equation to the equations set 
        /// </summary>
        /// <param name="eqv"></param>
        public void addEquation(double[] eqv)
        {
            if (_numOfVars == -1)
                _numOfVars = eqv.Length;
            else if (_numOfVars != eqv.Length)
                throw new System.IO.IOException("Inconsisstent number of variables in vector");

            _matrix.Add(eqv);
        }

        /// <summary>
        /// Clears the equations set
        /// </summary>
        public void clear()
        {
            _numOfVars = -1;
            _matrix.Clear();
        }

        private void sortMatrix()
        {
            for (int i = 0; i < _matrix.Count; i++)
            {
                if (_matrix[i][i] == 0)
                {

                    for (int j = i + 1; j < _matrix.Count; j++)
                    {
                        if (_matrix[j][i] != 0)
                        {
                            // swap
                            double[] temp = _matrix[j];
                            _matrix[j] = _matrix[i];
                            _matrix[i] = temp;
                        }
                    }      
                }
            }
        }

        /// <summary>
        /// Normalizes the vector according to the element on the specified position
        /// </summary>
        /// <param name="vector">Vector to be normalized</param>
        /// <param name="position">Position of element of normalization</param>
        /// <returns>Normalized vector</returns>
        private double[] normalizeVector(double[] vector, int position)
        {
            double[] result = new double[vector.Length];
            double factor = vector[position];
            if (factor == 0) return vector;

            //Normalize
            for (int i = 0; i < result.Length; ++i)
                result[i] = vector[i] / factor;

            return result;
        }

        /// <summary>
        /// Solve the equation set
        /// </summary>
        /// <returns>Solution vector</returns>
        public double[] solve() 
        {
            double[][] matrix = new double[_matrix.Count][];

            sortMatrix();

            //Copy list to array of arrays
            for (int i = 0; i < _matrix.Count; ++i)
                matrix[i] = _matrix[i];


            //Scan each row
            for (int row = 0; row < matrix.Length; ++row)
            {
                sortMatrix();


                //Normalize
                matrix[row] = normalizeVector(matrix[row], row);

                //Scan each column
                for (int otherRow = 0/*row + 1*/; otherRow < matrix.Length; ++otherRow)
                {
                    if (otherRow == row)
                        continue;

                    if (matrix[row][row] == 0)
                    {
                        //throw new Exception("0");
                    }
                    double multiplier = matrix[otherRow][row] / matrix[row][row];
                    
                    //Scan each cell
                    for (int column = row; column < matrix[row].Length; ++column)
                    {
                        matrix[otherRow][column] = matrix[otherRow][column] - multiplier * matrix[row][column];
                    }
                }
            }


            //Prepare the solution vector
            double[] result = new double[matrix.Length];
            int length = result.Length;
            List<double[]> resultMatrix = new List<double[]>();

            for (int i = 0; i < length; ++i)
                result[i] = matrix[i][length];

            return result;
        }
    }
}


